Abstract. A microscopic, quantum mechanical model for neutron vortices in the crust 
of a neutron star is presented. After a brief introduction to the Bogoliubov-de Gennes 
equations, which form the basis for our calculations, we present results for density distri- 
butions, vortex core sizes and vortex energies, both for an isolated neutron vortex and for 
the case when the vortex core overlaps with a cylindrical nucleus. Earlier results on the 
vortex core size are confirmed, indicating a much less dramatic variation of the vortex 
core size with density than predicted by the BCS formula. 
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1. Introduction 

A natural probe of the internal stucture and dynamics of neutron star crusts is rep- 
resented by pulsar glitches, which arc sudden accelerations of the star's rotational fre- 
quency. The observation of more than seventy events in about thirty pulsars provides 
a general view of the phenomenon. At present the standard model to explain pulsar 
glitches focuses on the interaction of superfluid vortices with the crust of the neutron 



star (Pines & Alpar(1985)). Some variants include crust breaking, resulting in a star- 
quake, or stellar plate tectonics due to the action of magnetic threading of proton vortex 
lines. 

The standard model is built in analogy to type II superconductors, where the attrac- 
tive interaction of vortex cores with the lattice defects of the metal prevents the motion 
of vortex lines through a metal sample. An important but unknown parameter is the 
pinning energy of a vortex to a nucleus, the difference in energy between the configura- 
tions in which a vortex and a nucleus are superimposed and the one in which they are 
far apart and non-interacting. Usually this quantity is calculated simply as the loss of 
condensation energy between the two configurations. Since the core of a vortex is com- 
posed of normal matter and the neutron superfluid which occupies the nuclear volume 
is, for sufficiently high densities, a weaker superfluid (i.e. has a smaller local pairing gap) 
than the superfluid outside the nuclei, it is energetically favourable for a vortex to remain 
attached to a nucleus. However, this analysis neglects important effects like the change 
in kinetic energy (Epstein & Baym 1988; Pizzochero, Viverit & Broglia 1997) and the 



presence of bound states (De Blasio & Elgar0y(1999)) 



Furthermore, the vortex is depicted as a cylinder composed of uniform normal mat- 
ter with a sharp radius equal to the BCS coherence length of the superfluid (the BCS 



coherence length is the mean square radius of Cooper pairs, see De Blasio et al.(1997) 
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for the case of hadronic matter). This is a very rough model since the gap changes 
smoothly inside the vortex and not abruptly, as assumed. The vortex-nucleus interaction 
results from effects requiring a full quantum mechanical solution of the BCS equations 
in the cell. Another approach which has been attempted is the Ginzburg-Landau model 
( Epstein fc Baym(1988)[ ). This model provides a good treatment of the kinetic energy 
density, which is an important ingredient in the energy balance within the interacting 
volume. Unfortunately, the pairing is not well treated in the Ginzburg-Landau equations 
since the two main conditions for their applicability, that the temperature should be close 
to the critical one and that there should be only slow and smooth spatial variations in 
the order parameter, are not fulfilled in a neutron star crust. 

The relevant quantity, the vortex-nucleus pinning force, can in principle be calculated 
from the variation of the pinning energy with respect to the vortex-nucleus distance. An 
approximate expression can be found by dividing the pinning energy by the minimal 
distance where the two objects can be considered as separate. The pinning force enters 
in two important ways in the standard model for pulsar glitches. Vortices are released 
from the crust when the frequency lag between the superfluid and the crust reaches a 
critical value. They unpin in clusters, decreasing suddenly the angular momentum of the 
superfluid. The glitch results from the conservation of angular momentum, causing the 
crust to increase its spinning rate. Therefore, the critical frequency is proportional to the 
value of the local pinning force. 

The second issue where pinning force is relevant is postglitch relaxation. This is 
particularly interesting, because time analysis following a glitch is now available for some 
glitches and especially for the Vela pulsar. In the vortex creep model, vortices that have 
abandoned their pinning centers repin with other nuclei within characteristic times which 
also depend on the value of the pinning force. The pinning force is of the order -Ep/^, 
where ^ is an appropriate length scale of the vortex-^nucleus interaction potential, roughly 
the size of the vortex core, and Ep is the pinning energy. This length is usually taken to 



be the BCS coherence length, but recent calculations ( De Blasio fc Elgar0y(1999) ) show 
that the vortex core can be significantly smaller than the BCS coherence length. This 
indicates that if the pinning energies were comparable, the pinning force should be larger. 
It is clear that a better understanding of the pinning model for pulsar glitches will be 
possible when more reliable calculations of the structure of a superfluid neutron vortex 
are available. 

Again, the study of superfluid vortices in neutron stars is not exhausted by considering 
solely neutron vortices in the crust. Below the inner crust, neutrons couple in the 
channel to generate a triplet superfluid, while protons, having a much smaller density, 
couple in the singlet channel. The flux of protons around the vortex core generates a 
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magnetic field parallel to the rotational axis. In addition, a large number (10^^ cm~'^) of 
proton flux lines is generated by the very high magnetic field (~ lO^^G) inside the core. 
The interaction of these flux lines with neutron vortices drifting towards the exterior 
of the star might be an efficient mechanism for magnetic field decay (see for example 



Ruderman(I997)[ ). 

The study of these physical situations may require a more accurate model for vortices 
in hadronic matter. A fully quantum mechanical model of vortices in type II supercon- 



ductors has been extensively investigated in some publications (|Gygi fc Schliiter(199I)| , 



Hayashi et al.(1998) ). It is possible to study the properties of neutron vortices paral- 
lelling the model for superconductors. Experiments with scanning-tunnelling microscopy 



and refined numerical calculations ( Gygi fc Schliiter(199I) ) have confirmed the theo- 
retical prediction Caroli, de Gennes fc Matricon(1964) , Bardeen et al.(I969) that bound 
states are formed in the center of the core of a vortex in type-II superconductors. The 
pairing gap decreases to zero at the center of the core and grows to an asymptotic value 
within distances of the order of or larger than the coherence length. 



In the present paper we extend our previous work (De Blasio & Elgar0y(1999)), here- 
after Paper I, on the microscopic structure of neutron vortex lines. In Paper I we made 



use of the Bogoliubov-de Gennes equations (de Gennes(I989), Ketterson & Song(1999)), 
that have been successfully developed in studies of vortices in type II superconductors 



( Gygi fc Schhiter(I99I) , Hayashi et al.(1998) ) and more in general non- homogeneous su- 
perconductivity. 

In Paper I we were mostly interested in the vortex core size, since this quantity influ- 
ences the pinning energy and force. We will here continue our study of this quantity, pro- 
viding more complete results and extending our calculations towards a fully microscopic 
calculation of pinning energies. Furthermore, we will show results for quantities such as 
density distributions, pairing potentials and vortex tension. The calculations presented 
here are also more refined than those in Paper I, as they are fully three-dimensional 
and include effects like the Hartree mean field and a density-dependent pairing force. 
Furthermore, since our earlier paper was brief on the technicalities, we will give a more 
detailed account of the formalism and the numerical solution of the relevant equations. 



2. The model 

The ^Sq neutron superfluid in the inner crust is spatially nonhomogeneous due to the 
presence of a nuclear lattice and the vortices induced in the superfluid by the rotation of 
the star. To describe this system one needs to go beyond the standard BCS formalism 
for fermion pairing. A formalism for this has existed in solid state physics for several 
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years in the form of the Bogohubov-de Gennes (BdG) equations ( |de Gcnnes(1989) ). For 
readers not famihar with these, we will give a sketch of their derivation. 

The theory is most conveniently formulated in terms of the neutron field operators 
and ipti^) which respectively destroy and create a neutron with spin projection a 
in position r. As a starting point, we choose a pure pairing Hamiltonian 
H = Ho + i/int = Ho- \g\ J dV^jW^jW^iWV'TW, (1) 
where \g\ is the pairing strength (we use the notation \g\ to avoid confusion about the 
sign of the pairing strength, which in some works is defined to be negative, in others to 
be positive). This corresponds to a zero-range pairing interaction of the form 
i;(r,r') = -|5l'5(r-r'). (2) 
and necessitates the use of a momentum cutoff confining the interaction to a narrow 
set of states near the Fermi surface. Note that in general the pairing strength is space 
dependent due to density variations induced by the mean field or by the vortex. We will 
return to this point later, for the time being we assume that it is a constant. 

The one-body hamiltonian has the form 
Ho^ J d\[i;]{r)hoi;^{r)+i;l{v)ho^Pi{r)] (3) 
where 

ha^i+W{r) + S{r)=i+Y{r) (4) 
and 

2 m 

In the above equations A is the chemical potential, W^(r) is an external potential repre- 
senting a nucleus, S{r) is a self-consistent potential generated by the neutrons themselves 
while Y{r) stands for the total field. 

By selective averaging, or, more formally, through applying the saddle-point approx- 
imation to the grand canonical partition function which follows from this Hamiltonian, 
one can derive a mean field Hamiltonian 

73^L/,t/'„\,/,t/' 



Hmfa = Ho + J dMi'l (r)V'I (r) A(r) + A* (r) (r) Vt (r)] + C (6) 
with 

A(r) = |g|(^i(r)^TW> (7) 
being sometimes called "the complex order parameter" , "the pairing potential" , or "the 
pair function" . Since this is a mean-field approximation, quantum fluctuations in the 
order parameter is neglected. The quantity 

C=/.V^ (8) 

is a c-number which plays no role in the free energy minimization leading to the 
Bogoliubov-de Gennes equations, but which is important in the energy calculations. Here, 
and in the following (. . .) denotes thermal averages. The mean field Hamiltonian (ph can 
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be diagonalized by a canonical Bogoliubov transformation 

a,,T it) = J d'r[U* (r) (r , i) - V* (r) V| (r, t)] (9) 

<iW = J d-'r[U,{ry^\ir,t)+V,{r)i;^{r,t)], (10) 
with the inverse transformation being 

Mr,t) = Y.iU^{r)a,^^{t) + V:{r)alJt)] (11) 

i 

V|(r,t) = Y,iU:{M,iit) - Viir)a,,^{t)]. (12) 

i 

In these equations, i is a label for the single-particle eigenstates. The coefficients U and 
V of this Bogoliubov transformation must obey the canonical relations 
5^[C/,(r)[/;(r') + Vi{r')V;{r)] = 6{r - r') (13) 

i 

J2mr)V;{r')-U,{r')V:{r)] = (14) 

i 

I cfr[Ui{v)U*,{v) + y,(r)y;(r)] = Sw (15) 

d''r[U,{r)VAr) -UArmir)] = 0, (16) 
in order for the quasiparticlc operators to follow Fermi statistics. Carrying out this trans- 
formation, the Hamiltonian takes the simple form 

i 

The constant term Uq is the ground state energy of the system. The Bogoliubov-de 
Gennes equations for the transformation coefficients which bring the Hamiltonian to the 
diagonal form Hb follow from the Heisenberg equations of motion for the field operators, 

i ip= [^jiJg], that is, in matrix form 

,d_ I V'T(r,t)\ ^ f i+Y{r) -A(r) \ / i^^{r,t)\ 

The i/i-operators are not eigenvectors of the Hamiltonian and thus have no definite fre- 
quency. We can, though, express them through the eigenoperators of the Hamiltonian, 
a, a^: 

ai,T(t) = a,,te-*^'*; 4^(t) = al/^"'. (19) 
By doing so, and gathering the terms with aq^-^ and so on, one obtains the BdG equations 

a-(r) -i-Y{,) ]\Vt) \Vt ) 

Solving these equations, one gets both the amplitudes U and V along with the excitation 

energies Eq of the system. Once the amplitudes are known, the expectation value of any 

single-particle operator O can be obtained from 

(O) = 2 ^ / <fr[U;OUim) + ViOV:{\ - m))], (21) 

Ei>0-' 

where f{Ei) = {l + ex.p{Ei/kBT)~^ is the Fermi distribution function, fcs is Boltzmann's 
constant and T is the temperature. More explicitly, we write the BdG equations as 

"1^ + ^(r) - Ef) Ui{r) + A{r)Vi{r) = EiUi{r) (22) 
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, y(r) - Ef mir) + A*(r){/,(r) = E,V^{r) (23) 
2m J 

where U and V are the quasiparticle ampUtudes, Ep is the Fermi energy, A(r) is the 
(space-dependent) pairing potential and m is the neutron mass. In the above equations 
we have replaced the chemical potential with the Fermi energy. The resulting violation in 
particle number conservation is found to be very small. In equations (^ 23 ) the subscript 
i represents all relevant quantum numbers. The pairing potential has to be calculated 
self-consistently as 

A(r) = |g| J2 U.ir)V:ir){l-2f{E,)) (24) 

i;0<\Ei\<hn 

where the sum is over quasiparticle states having energy eigvenvalues with absolute values 
\Ei\ smaller than a cutoff hfl, in our calculations taken to be 10 MeV. In general the 
set of quantum numbers i includes the principal quantum number, the projection of the 
angular momentum /i — 3/2,— 1/2,1/2,3/2,... and the wavenumber fc^ parallel to 

the z-axis. As we want to describe neutron vortices with cylindrical symmetry, we write 
the quasiparticle states as 

Ui{r) = ^u„^fe^(p)exp[i(/x- 1/2)0] exp[i(fc^z)] (25) 
VL 

Vi{r) = -^w„^feJp)exp[i(/^ + 1/2)6*] exp[i(fc^z)] (26) 
where p, 9, z are the cylindrical coordinates, L is the length of the cylinder, n is a radial 
quantum number and [i is half an odd integer. 

The contributions to the pairing potential are calculated in each subspace of fixed /i 
and kz and are then summed to give the total pairing potential. In Paper I we considered 
the simplified case of a pairing potential calculated at fcz = 0. This approximation is the 
one currently used in self-consistent calculations of vortices in type-II superconductors, 
where the electron effective mass is anisotropic and increases strongly along directions 



perpendicular to the vortex axis ( Gygi fc Schhiter(1991) ). In the case of vortices in neu- 
tron matter, variations in the effective mass are expected to be much smaller. In the 
present study we shall therefore go beyond this approximation and allow for different 
values of the wavenumber along the vortex axis. Since a direct comparison of vortices in 
neutron matter and in superconductors is important, we shall also refer to the case with 
fcz = whenever such a comparison can be revealing. 

The angular dependence in equations ( p5| , p6| ) of the quasiparticle states follows from 
imposing a pairing gap of the form 

A(r) = A(p) exp[-i0] (27) 
representing a vortex with one quantum of circulation. We are not interested in 
larger winding numbers, since they are expected at much higher energies. Following 
( |Gygi fc Schliitcr(1991) ) we expand the quasiparticle states in terms of cylindrically sym- 
metric Bessel functions and impose the boundary condition A(i?) = at the edges of a 



= (31) 
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cylinder of radius R, where R is typically 50-100 fm. More specifically, the basis functions 
are chosen as 

<^jm(p) = n J ^ \Jm ("5^4) ' j = 1, ■ ■ • , A^, (28) 

where m = /i± ^ is an integer and aj„i is the jth zero of the Bessel function J,n{x). The 
dimension N of the basis is chosen large enough to ensure convergence and stability of 
the quantities of interest. The quasiparticle amplitudes are expanded as {i — {n/ifc^}) 

Ut{p) = ^Qj0^^_i(p) (29) 

3 

Vtip) = Yd^j(t)j^+^{p). (30) 

3 

We consider first the case without mean field, which is the one usually addressed in the 
calculations for type-II superconductors. For given values of /i and equations ( |2^p3| ) 
can be written as a 2N x 2N matrix eigenvalue problem 
T- A 

AT _y+ 

where the superscript T denotes the transpose of a matrix, = 

(c„i, . . . , C„jVj dnl, ■ ■ ■ , dnN), 

and the matrix A is given by 

Aji' = / 0jM-i/2(p)A(p)0y^+i/2(p)pdp. (33) 
Jo 

In cases where a mean field is present the matrix elements in equation ( |3^ ) get an 

additional term. In our calculations we assume that the mean field of a nucleus has a 

Woods-Saxon shape, 
Wn 

l + exp(^) 

where typical values of the constants are Wq ~ —30 to ^ —40 MeV, Rn ^ 5 fm, and 
a ~ 0.6 fm. Note that we only consider cylindrical nuclei in the present paper. The case 
of a spherical nucleus interacting with a cylindrical vortex line is quite complicated since 
the two objects have different symmetries. However, it is worth noting that several cal- 



culations (Lorenz, Ravenhall & Pethick(1993), Pethick & Ravenhall(1995) ) predict that 
nuclei take on cylindrical shapes in parts of the inner crust of a neutron star, so our 
results have in fact a direct bearing on the actual problem. 

Even without a nucleus present, the neutrons move in a self-consistent Hartree field 
of the form 

Sip) ^ -\g\nip) (35) 

where n{p) is the neutron density distribution within the cylinder, given by 

n{p) = 2L-' J2 \vn.,.kAp)\'- (36) 
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In the presence of a mean field with cyhndrical symmetry, Y{p) = W{p) + S{p), the 
matrix elements of the one-body operator in equation (^2|) acquire an additional term of 
the form 

W%,=T rfw<^j^±i/2(p)y(p)'/'j>±i/2(p). (37) 
Jo 

For generality the expressions have been written at finite temperature, making use of the 
Fermi occupation factors for quasiparticle excitations. Although the case at finite tem- 
perature presents no additional difficulty, in applications to neutron star superfluids the 
temperature is essentially zero compared to the Fermi temperature. Thus, all numerical 
calculations in this paper are carried out at T = 0. One can use the symmetries of the 
BdG equations under the transformation 

{U,V,E)^{-V*,U\-E), (38) 
i.e., if {U, V, E) is a solution, then so is {~V* , U* , — i?), to rewrite equation ( p^ simply 
as 

A(r) = 2|g| ^ U,{v)V* (v) = A'^' ^ u^.^^kMvn.^.kJp) (39) 

i;0<Ei<hn n,fj,,k^;Q<En^fi^k^<fi^ 

where the sum now extends over positive eigenstates of the BdG equations only. This 
symmetry can be in principle exploited to solve the system (^l]) for only positive (or only 
negative) eigenstates. 

A point-like pairing interaction v{r, r') = — |g|(5(r — r') has been used in previous work 



on superconductors (de Gcnnes(1989), Gygi & Schliiter(1991), Hayashi et al.(1998)). In 



the case of neutron matter this approximation should work whenever the range of the 
inter-particle interaction (which is of the order ^ 1 fm) is at least comparable to or 
smaller than the nucleus and the vortex. Exactly as for a system of electrons, a cutoff 
energy hfl has to be introduced for the sum in equation ( |39| ) to converge. 

3. Results 

3.1. Calculations for different gaps 

One of the parameters still largely unknown in the physics of neutron stars is the value 
of the neutron pairing gap as a function of the density in the crust and in the interior. 
The reason is that the gap is very sensitive to the value of the parameters defining the 
neutron-neutron interaction at the Fermi surface, as can be seen from the weak-coupling 
formula for the gap 

Aw2?il7exp( I (40) 

where hft is an energy cutoff (often taken to be the Fermi energy), iV(0) is the density 
of states for a fixed spin projection at the Fermi surface and V is the effective interac- 
tion between two Cooper-paired neutrons. In many studies (Baldo et al. 1990; Elgar0y 
& Hjorth- Jensen 1998; Khodel, Khodel & Clark 1996) the last quantity is approximated 
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by the free neutron-neutron interaction, taken from various nucleon interaction models 
fitted to proton-proton and neutron-proton scattering data. The weak-couphng approx- 
imation is poor in neutron matter due to the strong momentum-dependence of the free 
neutron-neutron interaction, and the BCS gap equation has to be solved in its full com- 



plexity (Baldo et al.(1990)). While sometimes numerically tricky, this task can now be 
carried out accurately, and the results turn out to be insensitive to the model cho- 
sen for the nucleon-nucleon interaction as long as it fits the scattering data accurately 



( [Khodel et al.(1996)i |Elgar0y fc Hjorth-Jcnsen(1998)| ). However, the interaction between 
two neutrons is modified when other neutrons are present, and an important contribution 
comes from the so-called induced interaction, corresponding to the exchange of density 



and spin-density fluctuations between two neutrons ( Migdal(f 967) ) . The effect of this has 
been calculated by Wambach, Ainsworth, & Pines(1993) and Schulze et al.(1996) with 



somewhat different results. Wambach et al.(1993), calculating the induced interaction 
accurately, but using the weak-coupling approximation for the gap, predict a decrease by 
a factor of ^ 3 at all densities compared with the calculations using free nucleon-nucleon 



interactions. Bchulzc et al.(1996), solving the full gap equation, but using a somewhat 
rougher treatment of the induced interaction, predict a similar reduction of the gap at low 
densities, but an increase at higher densities. We know of no estimates for contributions 
beyond the induced interaction; some of these may turn out to be important. Given the 
uncertainties in the predictions for gaps in uniform neutron matter, we will in some of 
our calculations let this be a free parameter and study how other quantities depend on 



the size of the gap. In other calculations we will use the results of [Wambach et al.(1993)| , 
since these are commonly employed in other neutron star studies, and have so far been 
consistent with the glitch data. Actually, we find the nice feature that some important 
properties of the vortex are not too sensitive to the gap. In a first group of calculations 
we shall thus keep the value of the gap free, using several values of the pairing force \g\ 
for a fixed Fermi wave number. 

Keeping comparable values of the gap at infinity and changing kp wc find that |f;| 

is approximately inversely proportional to the Fermi wave number. This is in rough 

agreement with the weak coupling formula which predicts (with V{kF, kp) — -~\g\) 

, , 271^^2 1 ^ ^ 

i.e., \g\kp « constant modulo a slowly varying logarithmic factor. 

Figure |^ shows the pairing gap as a function of the distance from the vortex core. As 
analyzed in Paper I, the gap increases from zero to an asymptotic value Aqo, which is 
the value in infinite homogeneous matter. From the variation of the gap we can extract 
the appropriate size of the vortex. As in Paper I, we take the intersection of the tangent 
line at the origin p — Q with the line Aqo to define the length ^2- The quantities ^50 
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and are the distances from the axis where the gap reaches 50% and 90% of the value 
at infinity, respectively. The length is simply the BCS coherence length, i.e. the size 
of a Cooper pair. In Paper I we showed that with one single = state the size of 
a vortex tends to change smoothly as a function of the pairing force \g\ and the Fermi 
wave number kp. Table 1 of the present paper collects data for the vortex size for several 
different Fermi wave numbers and asymptotic gaps. The data show that the vortex size, 
although certainly not constant, has a limited range of variation, while the BCS coherence 
length changes dramatically. For a fixed kp and increasing Aqo the lengths generally tend 
to decrease. The reason for the different scaling of the vortex core compared with the 
BCS coherence length will be examined in Sect. 3.2 while consequences of astrophysical 
relevance will be discussed in Sect. 3.6. 

Figure ^ shows the eigenvalues of the vortex state for the case of one single kz 
mode. Similar to that found in other systems exhibiting fermionic superfluidity, such 
as superconductors ( Ketterson fc Song(1999) ) or finite nuclei ( ^chuckfc Taruishi(1996) ), 



a window of width 2Aoo opens up where only bound states are present. These are 
visible as a branch that for large angular momenta /i approaches the energy of con- 
tinuum states, very close to what one finds for vortices in type II superconductors 
( |Gygi fc Schliitcr(1991) ). The smallest eigenvalue has an energy of the order ~ A'^/Ep- 



For more k^ states, as shown in figure ^, there are in general more bound states for a 
given angular momentum ^ and the energy is roughly E ^ /iA^/_Ei?[l — (kz/kp)'^]^/'^ 



(Caroli, de Gcnnes fc Matricon(1964)). The eigenfunctions of the bound states are con- 
fined in the region of the vortex core and decrease exponentially far from the core. The 
wavefunctions of the continuum oscillate like scattering states with a radial wave number 
kp determined by the condition E^^ ~ + {h? /2m)'^{k'j, + fc^ — fc|,)^. 

3.2. Calculations with fixed gap at infinity 

In a second group of calculations we assume that the pairing gap of neutron matter 



is given by the results of (Wambach et al.(1993)). In short, once the neutron Fermi 



wavenumber is fixed, we choose the value of the pairing strength g which reproduces 
the given value of the pairing gap. This procedure is necessary when a nucleus is present 
in the center of the vortex core, because a nucleus modifies the local density. This in turn 
has strong influence on the local gap and condensation energy density, which are both 
strongly density-dependent. 

We solve the BdG equations for a homogenous system at various densities, and require 
that the energy gaps thus obtained be equal to the ones of Wambach et al. This provides 
us with a density dependent pairing strength \g\ — \g[n{p))\. When solving the BdG 
equations for a vortex line, we naturally find that the neutron density varies within 
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the cylinder. Thus, g will also vary throughout the vortex, and this must be taken into 
account. The pairing potential is then modified to 

A(r) = 2 l5(r)|t/.(r)F;(r), (42) 

v,o<Ei<nn 

and similarly for other quantities involving |f;|. 

In the numerical solution of the BdG equations, we use 50 — 100 Gaussian mesh 
points for the radial coordinate, and some 50 — 100 angular momentum states and 8 — 10 
plane wave states (for the z direction) in the expansion of the U and V amplitudes. 
Starting from initial approximations to A(/9) and n{p) the BdG equations are solved 
by diagonalizing the resulting eigenvalue problem. This gives us approximations to U 
and V , and from these new approximations to A(/9) and n{p) are obtained. We iterate 
this procedure until the variations in these quantities are small from one iteration to 
the next. Usually this procedure converges after 5-10 iterations, and it is adequate to 
use A(/9) — constant and n{p) — constant as initial approximations. Figures I and I 
show the pairing potential A(p) and the neutron number density n{p) respectively as 
functions of the distance from the axis for kp — 0.8 fm^^ in four different configurations: 
an isolated vortex at p = 0, an isolated nucleus, vortex and nucleus both present, and 
uniform matter. When only the vortex is present the gap rises linearly to an asymptotic 
value within lengths specified by the parameters ^ defined in Section 3.1. The fact that 
the pairing strength is now space-dependent due to the decrease of the density in the 
core visible in figure ^, changes only slightly the values of the vortex core radius. The 
shape of the gap is also only slightly altered by the variation of the pairing strength. 
The case where only the mean field is present shows a decrease of the pairing gap in the 
region occupied by the nucleus, qualitatively consistent with semiclassical approximations 



(Broglia et al.(1994)). When both a nucleus and vortex are present the picture becomes 
more complex due to the strong distortion of the velocity field induced by the nucleus. 
To understand the difference in behavior of the gap in the presence of a vortex, with 
or without the nuclear mean field, we analyze the problem using the Ginzburg-Landau 
(GL) equation 

- — 7/;(r) + VF(r)V'(r) - AV^(r) - B|V'(r)|V(r) = (43) 

4m 

where ^{v) is the order parameter and the self-consistent mean field S{p) is neglected 
compared to W{p)^ which is a very reasonable approximation. The parameter A is nega- 
tive for T < Tq, where Tq is the critical temperature of the superfluid, while B is always 
positive. Using cylindrical coordinates, taking i/j{r) = e^*^/(p), W{r) — W{p), and using 
units where /Am = 1, equation (E3|) can be written as 



np) + -/{p) 



-W{p)~A~B\f{pt-^^ 



fip) = 0. (44) 



Elgar0y and De Blasio: Superfluid vortices in neutron stars 



13 



The boundary condition at p = is /(O) 0. Looking at small p, we neglect the non- 
linear term, and distinguish between two cases. First, for the case of no nuclear mean 



|) becomes 



field, equation 

f'ip) + -fip) ^f{p) = 0, (45) 
P P 

since A is supposed to be small near Tc and thus negligible compared to for small 
p. This equation is easily seen to have the solution f{p) cx p. Thus, we expect the pairing 
potential to show a linear increase for small p in the case of an isolated vortex, and this 
behaviour is also seen in our microscopic calculations. Looking at the second case, an 
attractive mean field W{p) ~ —Wq^ assumed to be constant for the limited range in p 
we are looking at now, equation (^ is reduced to 

f'ip) + -/(P) + (wo - ^) fip) = 0, (46) 
where we have absorbed the constant A in the mean field Wq. Introducing a new variable 



through p ■ 



/Wo 



we obtain 



x^f'ix) + xf{x) + {x' - \)!{x) = 0, (47) 
which the reader may recognize as the Bessel equation of order 1. The solution satisfying 
/(O) = is f{x) oc Ji{x). Thus, in the case of a vortex interacting with a nucleus, we 
may expect the pairing potential to show oscillations at small p. 

Figure ^ shows the relevant lengths for the case of polarized gap. The lengths remain 
limited (< 15 fm) for essentially all densities and tend to decrease at higher density. 
This behavior is different from the BCS coherence length which scales like kp/ A^. In 
the realistic case where the pairing strength is fitted to the gap in uniform matter, the 
coherence length ^2 decreases as a function of fc^- Thus, the density dependence of the 
pairing strength \g\ is also very important, and the vortex size appears to depend in a 
quite complex way on both \g\ and the Fermi wavenumber. In the figure we also show the 



analytic estimate of the vortex core size made by Kramer & Pesch(1974) whose value in 



our units (energies in MeV and lengths in fm) is 
_ 661.3 



(48) 



Since the pinning force between a nucleus and a vortex behaves like ~ 
an increase of the pinning force of the order ^ Co/$2 can be expected. In addition we 
can confirm what we found in Paper I, that since the diameter of the vortex is always 
smaller than the lattice spacing, the vortex envelops at most one nucleus along a plane 
perpendicular to its axis. 

Another interesting parameter is represented by the coherent flux of neutrons in the 
presence of a vortex. The flux is given by 



J(r) 
Jip) 



h 

2mi 



E 



y,(r)VF;(r)-t/;(r)Vy,(r) 



which in the present case becomes 

^ ^^(/i-l/2)(|^.,(p)p- ^ 

E,<0 Ei>0 



2mp 



^{p+l/2)\v,{p)\' 



(7,(r)V(7*(r)-C/*(r)V[/.(r) 
(50) 



(l-/(i?,)) (49) 
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Figure |^ shows the flux of neutrons with and without a nucleus present in the middle 
of the cell for the case of kp — 0.4 fm^^. The main effect of the nucleus is to decrease 
the flux. The coherent velocity field can be obtained dividing the flux by the neutron 
density. In this case the velocity is even more depleted in the centre due to the much 
higher density in the nuclear region than outside, a situation typical for low values of kp. 

3. 3. Variations in A due to a center of scattering 

Pinning or antipinning between a vortex and a nucleus is due to local variations in 
both the kinetic energy density and the condensation energy. Although pinning will be 
examined in detail in a next section, it is interesting to study how the pairing gap changes 
in the cell when a nucleus alone is present. Some curves have already been presented in 
Figure 1, but here we shall be more systematic discussing of the applications to pinning 
calculations. The variation of the pairing gap as a function of the distance from the 
axis of a cylindrical nucleus when the vortex is absent can be seen in Figure |^ for some 
selected densities. Note that the presence of a mean field decreases the value of the gap. 
This is mainly due to the increase of the density in the nuclear region, since the pairing 
force l^l depends on the local density. Note also that the variation of the gap extends 
beyond the range of the mean field, indicated with a dashed line. This is due in part 
to the long tail of the density distribution beyond the radius of the potential, an effect 
due to the presence of scattering states. Secondly, there is a proximity effect between the 
nuclear region and the free neutron gas due to Andreev scattering (that is, scattering 
due to spatial variations the pairing potential) of quasi-particles. 



3.4^. Pinning energy and vortex self-energy 

The vortex self-energy or tension T, which is the energy carried by the vortex per unit 
length, is a relevant parameter for vortex-nucleus pinning. Let us consider the behavior 
of a single vortex line in interaction with the whole lattice. If the energy necessary for the 
vortex to reach the pinning centers is very high compared to the energy gained by pinning, 
the vortex responds stiffly to local deformations. The importance of this parameter can be 
grasped in a limiting case: for T 00 the vortex behaves like a rigid cylinder. For a vortex 
moving perpendicularly to its axis the total pinning force with the lattice arises from 
stochastic summation from all the pinning centers and is proportional to the square root 



of the vortex length (see for example Tinkham(1996) for the case of superconductors). 
Thus, the pinning force cannot balance the Magnus force, which increases linearly with 
the vortex length. If this was the case, the pinning mechanism would not be effective for 



storing the energy released during a glitch ( [Anderson fc Itoh(1975)[ ). For a finite stiffness, 
the situation is more complicated. Solving the equation of motion for a vortex line passing 
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through many centers of pinning would give the value of the total pinning energy which is 
able to be stored, but this is a difhcult task due to the difficulty in finding the appropriate 
boundary conditions in the presence of many centers of pinning. It has been argued that 
although finite, the tension is too large to account for large glitches ( Jones(1998)| ). If, on 
the other hand, the pinning model is correct, the value of the vortex tension is relevant 



for the dynamics of unpinning (Link & Epstcin(1995)) and for the value of the effective 
pinning force per unit length. 

The vortex tension is the energy difference per unit length between the configuration 
with and without a vortex, T — [i?(vortex) — i?(uniform matter)] /L. To calculate the 
energy of the superfluid in each configuration one can take averages of Eq. (6). The result 
is 

E= f d\^-^^ + 2Y, En-2Y,En f d'r\Vn{r)\\ (51) 

Table (2) shows the calculated points as a function of the Fermi wave number compared to 
the classical calculation where the kinetic energy density is integrated from the coherence 
length ^ up to an upper limit Rq: 



T = 4.38fc^ln[^yJ. (52) 
Evidently the presence of bound states does not change the vortex tension radically. 
Furthermore, constraining the vortex in a small cylinder can overestimate the role of 
bound states in the calculation of self-energy, because the vortex is macroscopic along 
the p direction. We also expect that uncertainties in the pairing gap or in the cutoff length 
i?o (that should be of the order of the vortex-vortex separation) may be much greater 
than the discrepancy we find between our microscopic calculations and the estimate 



(p2|). We conclude that equation (52) is reliable when applied to vortex calculations in 
a neutron star. Assuming a vortex-vortex separation of the order i?o ^ 10^^ cm, this 
formula implies a tension of the order T ^ 10^ erg cm~^ at kp — 0.8 fm^^. At such values 
a vortex can deform only at length scales of the order of several hundreds of fm, nearly 20 
times larger than the lattice spacing, unless the pinning energy assumes unrealistically 
large values. How a vortex line with such small deformability can adapt its shape to gain 



a net pinning energy is still an unsolved problem ( Iones(1997) ). We show in the following 
that an answer may be provided by the high value of the elementary pinning force. 

To calculate the pinning energy we need to consider four different configurations: a 
vortex superimposed on a nucleus, an isolated nucleus, an isolated vortex, and finally 
uniform matter. For each configuration, the energy needs to be calculated with Eq. (|5l|). 
The pinning energy per unit length can be calculated as 



Ep — 



/L. (53) 



£'(nucleus and vortex) — £^(nucleus) -I- £'(uniform matter) — £'(vortex) 
Table (3) collects some pinning data calculated with Eq. (49) and (51). A negative 
sign corresponds to a pinning situation in which it is energetically favorable for a vortex 
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to remain attached to a nucleus, while a positive sign indicates an antipinning situa- 
tion, where vortices are repelled by nuclei. We find an antipinning-pinning transition at 
kp = 0.30 — 0.40 fm~^. A similar transition was found in a Ginzburg-Landau approach 
( [Epstein fc Baym(1988)| . The pinning energy reaches values of the order of several MeV 
per unit length, corresponding to a very high value if compared to previous calculations 
performed solely on the basis of the condensation energy loss. As is also found in other 



models where the kinetic energy is accounted for ( Epstein fc Baym(1988) ), part of these 
large values can be attributed to kinetic energy effects. For such large values of the pin- 
ning energy a vortex is likely to deform sensitively in order to be able to catch the pinning 
centers. In addition, the ions in the lattice will shift to the vortex in response to the high 
pinning force, a situation referred to as superstrong pinning. Since nuclei in most of the 
crust are spherical, one might question the applicability of our calculations, where nuclei 
are considered to be cylindrical. Although the proper use of spherical geometry would, 
of course, quantitatively change the results, we do not expect dramatic differences. 



4. Conclusions 

The results presented in this paper could be a step towards an understanding of vortex 
structure and interaction in neutron star crusts. We have found differences between the 
predictions of microscopic theory and those of more macroscopic models. In particular, 
we showed that the BCS coherence length cannot be used as a reliable estimate of the 
vortex core size, especially not at low and intermediate densities. By solving the self- 
consistent problem of a vortex in neutron matter, we have shown that the vortex core 
has a much more complicated structure than that which is usually assumed in simplified 
treatments of the problem, e.g.: 1) the gap changes from zero at the axis of the core up 
to an asymptotic value, and 2): there are bound states present. 

We found that vortices in neutron star matter follow a behavior similar to that de- 



scribed by Kramer fc Pesch(1974) for superconductors, namely, at zero temperature the 
vortex core size scales as the inverse of the Fermi wave number. In a superconductor 
in clean limit this effect is associated with a linear increase of the vortex core size with 
temperature. We have not investigated the temperature dependence of our results, be- 
cause this is not very relevant in neutron star physics due to the low temperatures 
compared with the Fermi energy. Unfortunately, experimental data for clean type II su- 
perconductors are limited to the quite complex new cuprate systems. Scanning tunnelling 
microscopy has, however, confirmed the presence of bound states of quasiparticles in the 
vortex core, in good agreement with the Bogoljubov-de Gennes theory. In view of the 
small range of the nuclear force and of the simple structure of the Fermi surface, the 
BdG model could be a very good approximation for neutron star matter. 
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High pinning energies and small lengths imply a large pinning (or anti-pinning) force. 
In addition, the small value of the coherence length means that a vortex envelopes no 
more than one nucleus in a plane perpendicular to the vortex axis. The possibility for a 
vortex to store pinning energy is largely the result of a competition between the energy 
gain with pinning and the energy cost to deform. To be able to deform and catch a nucleus 
along its path, a vortex has to be relatively free to perform local deformations. Assuming 
that the energy of the vortex is solely due to kinetic effects, the tension turns out to be 
of the order T « 100fc|, MeV fm^i and is thus very high at densities corresponding to 
the deep inner crust. For a total pinning energy per nucleus Ep, a vortex can efficiently 
catch the pinning sites only if they are at least a distance « apart. For example, 

at kp = 1 fm^^ a vortex deformation of the order of the lattice spacing can occur every 
~ 8 or so lattice spacings, while at kp — 0.8 fm^^ the pinning efficiency would be 
much better. Calculations of the geometry of vortex deformation in a one-dimensional 



crystalline structure have been performed by (Link & Epstein(f 995) ). 



Finally, there are other unsolved problems in neutron star dynamics where a micro- 
scopic description of vortex states might prove helpful. The scattering of electrons off 
neutron vortices results in a coupling between protons and neutrons in the interior of the 
star. This is a possible coupling mode between the superfluid neutrons and the rest of 



the star, and depends on the microscopic structure of vortex lines ( Sauls (f989) ) 



The proton vortices in the interior of the star induced by the rotation will probably be 
similar to the neutron vortices examined in this paper, but additionally they will generate 
a magnetic field parallel to the rotational axis. How these magnetic vortices interact with 
the (much more numerous) flux tubes induced by the high magnetic field penetrating 
into the star is only one of the problems connected to the dynamics of superfluid vortices 
below the star's crust. Despite their importance in the evolution of the pulsar's magnetic 
field, the properties of proton flux tubes and their possible interaction with neutron 
vortices in the interior, where neutrons couple in a triplet phase, are not well understood 



( Rudcrman (1997) ). A microscopic description of proton flux lines might also be applied to 



the study uf ultra-high magnetic flelds and the evolution of magnetic field in magnetars. 
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Table 1. Gap at infinity Aqq and the length scales ^2, ^90 and 2^50 for different 
values of the Fermi wave number and pairing gaps. 





A 






S90 V^^J^i^j 




0.1 


0.0049 


344.9 


6.02 


11.11 


7.15 


0.1 


0.011 


153.4 


6.07 


19.89 


12.29 


0.2 


0.02 


152.0 


9.80 


15.61 


11.03 


0.2 


0.04 


81.83 


12.73 


16.56 


13.20 


0.2 


0.183 


18.41 


5.57 


8.04 


6.04 


0.2 


0.291 


11.57 


5.20 


7.13 


5.63 


0.2 


0.744 


4.54 


4.59 


6.18 


4.96 


0.5 


0.226 


37.34 


5.32 


17.71 


7.266 


0.5 


0.778 


10.71 


2.96 


9.13 


3.33 


0.5 


1.73 


4.87 


2.52 


6.75 


2.67 


0.5 


3.42 


2.46 


2.09 


3.24 


3.02 


0.8 


0.205 


65.8 


6.05 


24.60 


17.48 


0.8 


1.72 


7.84 


2.21 


9.21 


2.65 


0.8 


3.8 


3.55 


1.67 


5.02 


1.87 



Table 2. Vortex self-energy T (in MeV fm~^) at different Fermi wave numbers calculated 
from the program and the value Tf calculated from the analj^tical formula, where the 
coherence length ^2 has been used for the lower length b. The radius R of the cylinder is 
reported. 

kp (fm-^) T (MeV fm Tf (MeV fm"^) ^0 (fm) 

0.4 1.08 0.64 80 

0.5 1.03 1.344 80 

0.7 1.67 5.30 80 

0.85 3.06 7.401 60 
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Table 3. The pinning energy per unit length Ep (in MeV fm~^) at different Fermi wave 
numbers calculated for a potential well with the following parameters: radius Rn = 5 
fm, depth W = —30 MeV and diffusivity a = 0.6 fm. The radius R of the cylinder is 80 
fm. 



kp (fm-^) Ep (MeV fm-^) 



0.15 

0.20 

0.30 

0.5 

0.7 

0.85 

1.00 



0.68 

0.201 

-0.091 

-0.77 

-2.34 

-16.47 

-4.09 
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Fig. 1. Pairing potential A(p) as a function of the distance p from the symmetry axis of 
the vortex with a pairing force independent of the density and fixed to return the desired 
value of the gap at infinity. 
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Fig. 2. The eigenvalues of the BdG equations for a vortex in neutron matter for the 
case kp = 0.8 and a pairing gap at infinity Aoo=2.9 MeV. The positive and negative 
quasipartiele states in the continuum are separated by a gap equal to 2Aoo- The bound 
states are visible as the branch between the continuum states. Only one single value of 
fc = has been included in the calculation. 
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Fig. 3. Distribution of energy eigenvalues as a function of angular momentum, now also 
with different values for kz- Aqq = 0.5 MeV and kp = 0.22. 
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Fig. 5. Density distribution n{p) as a function of p for three different situations. The 
density of uniform matter is in all cases given hy kp = 0.8 fm~^. 
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Fig. 6. Various length scales for a vortex as a function of the density. The meaning 
the different symbols are explained in the text. 
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Fig. 7. The coherent velocity field for a vortex with and without a center of scattering. 
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Fig. 8. Pairing potential at various densities for the case of a nucleus immersed in a 
homogeneous (vortex-free) neutron superfluid. 



